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Summary 

Seismic waves are the most sensitive probe of the Earth's interior we have. With the dense data 
sets available in exploration, images of subsurface structures can be obtained through processes 
such as migration. Unfortunately, relating these surface recordings to actual Earth properties is 
non-trivial. Tomographic techniques use only a small amount of the information contained in 
the full seismogram and result in relatively low resolution images. Other methods use a larger 
amount of the seismogram but are based on either linearization of the problem, an expensive 
statistical search over a limited range of models, or both. We present the development of a new 
approach to full waveform inversion, i.e., inversion which uses the complete seismogram. This 
new method, which falls under the general category of inverse scattering, is based on a highly 
non-linear Fredholm integral equation relating the Earth structure to itself and to the recorded 
seismograms. An iterative solution to this equation is proposed. The resulting algorithm is 
numerically intensive but is deterministic, i.e., random searches of model space are not required 
and no misfit function is needed. Impressive numerical results in ID are shown for several test 
cases. 



1 Introduction 

Seismic waves provide the best and most direct probe available of the properties of the Earth's 
interior. In exploration seismology artificially generated waves recorded at the Earth's surface are 
commonly used to image material discontinuities in the subsurface (e.g. via migration). However, 
a more elusive goal is inversion, the determination of Earth properties such as wave velocity and 
density from seismic data. Tomographic inversions, which use only the traveltime information 
of isolated arrivals, are frequently used to obtain smooth, low resolution estimates of the wave 
velocity generally for use in imaging algorithms. Inversion methods which are designed to make 
use of all the information in the data, i.e., not only the arrival time and amplitude of pulses but 
the details of their shape as well, fall under the category of "waveform inversion" . The additional 
information contained in the waveforms promises to improve the accuracy and resolution of Earth 
structure models. 

Although elegant in theory, waveform inversion has not enjoyed much use partly because it is 
computationally expensive — more than an order of magnitude more expensive than migration — 
and partly because it requires expert user intervention to make it work. However, in more 
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recent years, faster computers, better inversion strategies (e.g. Pratt, 1999a, 1999b; Sirgue and 
Pratt, 2004; Yokota and Matsushima, 2004), and a greater demand for accurate information on 
subsurface structure have made this approach more and more appeahng. 

Most of the work done to date in waveform inversion falls into one or both of two broad 
categories: A) linearized methods and B) model-searching methods. Model-searching methods 
(see e.g. Sen and Stoffa, 1995) involve three main components: 1) the definition of a misfit 
or error function that measures the difference between seismic waveform data and synthetic 
seismograms for a proposed Earth model; 2) a pre-defined set of adjustable model parameters 
and parameter ranges; and 3) some rule or set of rules for searching the defined parameter 
space. The general idea is to find the combination of model parameters that minimizes the 
error. This sort of approach can work well if the number of possible parameter combinations to 
be examined is small. Otherwise, given that the generation of synthetic seismograms is usually 
quite expensive computationally in three spatial dimensions, searching a large-dimensional model 
space is impractical at this time. 

Linearized methods are, on the other hand, quite cheap computationally. A mainstay of this 
approach is to assume that the true, unknown Earth structure differs only slightly from some 
initial reference model. Under this assumption, one can assume that the true wavefield can be 
approximated by a combination of the wavefield in the reference model and a singly-scattered 
wavefield. This single-scattering assumption is known as the first-order Born approximation. 
Mathematically, this is written as 

V'(x, uj) ^ V^o(x, uj)+uj^ J dVG'o(x, uj; x')l^(x )^/'o(x', uj), (1) 

where ip and ipQ are the actual and reference wavefields, respectively; Gq is the Green's function 
of the reference medium; V is the difference between the reference and actual media (and is 
frequently known as the "scattering potential"); and u is the angular frequency. Because this 
approximation is linear in V, one can construct algorithms to extract V and, hence, an approxi- 
mation to the true model from measurements of on the Earth's surface. Of course, there are 
many practical considerations which can make this process difficult to perform, not the least 
of which is insufficient data. However, even from a theoretical perspective, there are problems. 
Specifically, if the true structure differs more than a little from the reference model — the usual 
case — the data can contain arrivals, most notably multiples, that cannot be predicted from the 
single scattering approximation. 

One approach to waveform inversion that does not fall into either the linearized or model- 
searching categories is what is known as inverse scattering. The name, "inverse scattering", 
comes from the knowledge that wavefields can be viewed as resulting from multiple scattering of 
wave energy within the Earth. The object is to find some way to invert these multiple scatterings 
to retrieve the Earth structure. The principal work done in this area uses an approach called the 
"inverse scattering series" . 

Originally developed by Jost and Kohn (1952) in quantum physics and later by Moses (1956), 
the properties of the inverse scattering series have been studied by Prosser (1980), among others. 
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Razavy (1975) was first to apply this method to the seismic inverse problem. In recent years, it 
has been most extensively explored by Weglein and collaborators (e.g. Weglein et al., 2003). The 
main idea of this approach is first to write the wavefield as measured on the Earth's surface as 
a Born series, the multiple-scattering generalization of eq. ([1]). In a compact symbolic notation, 
this is written as 

Ws = Ms + iGoV^o)s + iGoVGoV^o)s + . . • , (2) 

where the subscript indicates that we evaluate the wavefield only where we can actually measure 
it. In principle all orders of scattering off the unknown structure V are included in this sum. If 
we also assume that the quantity V can be written as an infinite series, 

V = Vi + V2 + Vs + ... (3) 

where is n-th order in the recorded data, we can insert this series into the Born series and 
collect terms of common order in the data. The result is an infinite set of equations for the Vn- 

= iGoV2^o)s + {GoViGoVMs 

: (4) 

These equations are meant to be solved in a cascade fashion, solving first for Vi then using that 
result to obtain V2 and so forth. 

By studying these equations numerically in one spatial dimension, Carvalho (1992) has found 
that this approach does not appear to be convergent for an arbitrary contrast between the 
reference and the actual medium. Following his result, Weglein and collaborators have studied 
the sub-series approach, in which various terms in the full inverse scattering series of common 
form are combined into separate sub-series, each of which is postulated to perform a different 
task corresponding to a conventional processing step (Weglein et al., 2003). Using the sub-series 
approach, they have been able to achieve some promising results with no obvious problems with 
divergence. 

We have derived an alternative approach to the problem of waveform inversion — one which 
also falls into the category of inverse scattering. We introduce this new development in the next 
section. 



2 Theory 
2.1 Derivation 

To introduce and test the concepts of our approach, we present the theory in one spatial dimen- 
sion. The data will consist of a single trace recorded at 2; = for a source also located at 2; = 0. 
The source-time function will be a 5-function, making this single trace the complete impulse 
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response for the chosen source/receiver pair. The unknown structure we will wish to recover is 
assumed to be located "beneath" the source/receiver, along the half-line z > 0. Further con- 
straints we impose on the physics of the problem are that we have purely acoustic propagation, 
fully variable velocity, no attenuation, and constant density. This last constraint is imposed not 
only for simplicity but for the reason that it is not possible to reconstruct both velocity and 
density from a single trace in ID. 

We take as our reference, or background, a constant- velocity medium with wave propagation 
velocity cq and associated Green's function Go- Let c{z) and G be the wave velocity and Green's 
function of the true medium, respectively. We define the scattering potential as 

y(^) - 4) - ^- 

The basis of our derivation will be what is known as the Lippmann-Schwinger equation, the 
integral equation equivalent to the more frequently seen differential equation for acoustic wave 
propagation: 

G{z,uj; zo) = Go{z,uj; zo) + ^ dz'Go{z,uj; z')V{z')G{z' ,uj; Zq). (6) 



One method of solving this equation results in the generation of the Born series or, at least, 
the generation of Born approximations of various orders. For instance, the first-order Born 
approximation (also known simply as the Born approximation) is obtained by approximating G 
by Go above on the right only, yielding 



G{z,uj;zo)^Go{z,uj;zo) + — dz'Go{z,uj; z')V{z')Go{z' ,uj; zq). (7) 

Cq J —oo 



The second-order Born approximation can be obtained by first replacing G on the right-hand 
side of eq. (P) with the entirety of the right-hand side. 



G{z,uj]Zo)=Go{z,uj]Zq) + — / dz'Go{z,uj]z')V{z')Go{z' ,uj]Zq) 

Cq J — oo 



/CO POO 
dz' dz"Go{z,uj-z')V{z')Go{z',u-z")V{z")G{z\u-zo), (8) 
y -oo J —oo 



and then setting G ~ Go on the right again. The general rules for obtaining approximations of 
arbitrary order are 



G^'^\z,u]Zq) = Go{z,uj;zo) 

G(")(^,^;^o) = Go{z,u;;zo) + ^ r' dz'Go{z,u;-z')V{z')G^^-'\z',u;zo), (9) 



roo 



C, 



where G*-"-* indicates the n-th approximant to G. For our purposes, however, the once-iterated 
equation without the substitution of Go for G, eq. ([H]), is what we want. 
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To greatly shorten some of the equations to follow, we define some helpful notation. Setting 
z = Zq = 0, we let 



(GoVGo)ii5 = — / dz'GoiO, u; z')V{z')Go{z' , u; 0) (10) 

and 

(GoVGoVG)fi5 = — dz' dz"Go{Q,uj-z')V{z')Go{z',u-z")V{z")G{z'\uj-Q), (11) 

Cn J —oo J —oo 



where the subscript RS indicates that the quantities are evaluated only for the given source/rec- 
eiver geometry. Thus, our once-iterated equations becomes 

G(0, u; 0) = Go(0, u; 0) + {GoVGo)rs + (GoVGoVG)^^. (12) 
Using the explicit form of the Green's function Gq, 

Go{z,u;zo) = '^e'^\^~^'\/^\ (13) 

and the fact that V{z) = for all 2 < 0, we find 

= -iv^V'(2w/c„), (14) 



where V is the Fourier transform of V. Inverting this relationship, we find 

4 /"CO 

TTCq J —oo 



4 r'^ 

Viz) = / duje-^'^'/^{Go\Go)Rs. (15) 

TTCn J —oo 



(Note that our Fourier transform convention is that of symmetric normalization for the inverse 
and forward transforms.) Letting D{uj) = (7(0, u; 0) — 6*0(0, 1^; 0) and 
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U{z) = / dLue-^''^'/''W{u;), (16) 

vrco J-00 



which is tantamount to a constant-velocity depth migration of the data, we apply the above 
inverse Fourier transform to the once-iterated equation and rearrange to get 



4 

V{z) = U{z) + — / ci^e-2-^/'=«(GoVGoVG)^s. (17) 

TTCq J -00 



Noting that G on the right-hand side depends on V , we see that this equation is an inhomoge- 
neous, highly non-linear Fredholm integral equation of the second kind for the scattering potential 
V . This equation is a novel result and forms the basis for our approach to waveform inversion. 
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2.2 An iterative method of solution 

There is no pre-existing mathematical apparatus for solving equations like eq. f|T7j) . One possible 
scheme for solving it — one that should eventually be vetted by mathematicians — can be obtained 
by analogy to the iterative way of solving the Lippmann-Schwinger equation that we discussed 
earlier (eqs. [H]). We propose to construct a sequence of approximants to V by the following rules: 

V^^\z) = (18) 
= U{z) + / dcje-^^^^/'^'MGoV^^-^^GoV^^-^^G^"-^)) , (19) 

where we generalize the notation of eq. (fTT!) by substituting ^("~^) for V and use G*^""^^ to 
indicate the Green's function for a medium with potential V'^'^~^\z). Going back to the once- 
iterated Lippmann-Schwinger equation, eq. (fT2|) . which holds for any medium, we have 

(GoV^^^GoV^^^G^-'))^^ = G^^HO, ^; 0) - Go(0, u- 0) - (GqV^-'^Go)^^ . (20) 

Using this relation in our iterative equation, we can achieve a substantial simplification to our 
second iteration rule above: 

y W [z) = U{z) + {z) + — r due-^''^'''' { ^("-1) (0, cu; 0) - Go(0, cu; 0)) . (21) 

In parallel with our definition of [/, eq. (fT6l) . we set 

(z) = / dcje-2*^"/^o ( G("-^)(0, cj; 0) - Go(0, cu; 0) ) (22) 

vrco J-oo ^ ' 

to get the final version of our iteration rules: 

V^''\z) = 

\/W(^) = U{z) + V^''-^\z)~U^''-^\z). (23) 

The algorithmic interpretation of these iteration rules is simple. At each step, the new 
approximant y*-"-* is obtained from the previous one, the "migrated" data U, and the migrated 
synthetic data from a medium with scattering potential V^("~^). Thus our algorithm should 
recover V through a deterministic sequence of forward modelling simulations. It is for this 
reason that we call our approach to inverse scattering "iterative inverse propagation" . 



3 Numerical Examples 



As proof of concept, we present three synthetic tests of ID iterative inverse propagation. In each 
case, we constructed a test model and generated synthetic seismograms with a staggered-grid 
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finite-difference (FD) code (see e.g. Virieux, 1986) using a Gaussian source wavelet with a half- 
width of about 0.5 s. In the inversion, we used an initial constant background velocity of 1 km/s, 
and the necessary forward modelling was done with the same FD code and source wavelet that 
was used to generate the synthetic data. We emphasize that no pre-processing, such as removal 
of multiples, was performed on the data. 
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Figure 1: One layer of thickness 3 km with a wave velocity 40% above background. Comparison 
of the true model (in red) and the inversion results (in black) are shown for various iterations 
indicated in each panel by the parameter A^. 

For numerical reasons, one deviation from the derived iteration rules, eq. (!23|) . was used. 
At each iteration, the model obtained from the previous step was used as the new "mapping 
velocity", i.e., the velocity used to map time into space. In other words, instead of using the 
factor of 2z/cq in the exponential of eq. fl23p . a more general function T{z) — the traveltime from 
the source to the point z and back to the receiver — was computed from the previous iteration. 
This very practical ad hoc feature was used to eliminate some instabilities that were otherwise 
occurring at the bottom end of the structure of any test attempted. In fact, this approach also 
accelerated convergence and indicates the need to incorporate non-constant reference models in 
future work. 
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Figure 2: A smoothed random model with velocity variations of approximately ±40% about 
background. Note that the inversion recovers even the gradient at the end of the model. 

Figure 1 shows the results for our first test, a model with a single high- velocity layer imbedded 
in an otherwise homogeneous medium. The results of the first three iterations and iterations 10, 
20, 30, and 40 are shown. Both the amplitude of the velocity jump (40% above reference) and the 
positions of the top and bottom of the layer are quite well recovered, with the essential features 
of the true model recovered by iteration 20 and only moderate improvement obtained thereafter. 
Those deviations from the true model that are present are due to the band-limited wavelet used 
in both the initial synthetic data and the FD modelling performed in the inversion. 

Figure 2 demonstrates the ability of the algorithm to reconstruct the details of a complicated, 
smooth model with fluctuations of approximately ±40% about the reference velocity. We see 
that by iteration 10 the errors in the inversion result are negligible. It is particularly notable 
that even the smooth gradient at the end of the model is recovered. 

Finally, Figure 3 shows the results for a single low-velocity layer of velocity 40% below ref- 
erence. The value of this example is that it shows clearly what happens to multiples as the 
inversion progresses. Specifically, one sees that they are gradually "pushed" upward into the 
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Figure 3: One low-velocity layer of thickness 3 km with a wave velocity 40% below background. 
Note the progressive collapse of multiples into the main structure. 



main structure, indicating that the information they provide is instrumental in obtaining an 
accurate final result. 



4 Discussion 



We have shown the potential of iterative inverse propagation to recover the details of wave velocity 
structures in ID. There are, of course, going to be many hurdles to overcome in extending this 
method to practical use. In addition to incorporating non-constant reference media, we must 
also accommodate irregular source/receiver geometries, elasticity, attenuation, limited frequency 
ranges, and unknown source wavelets, just to name a few. There will also be the issue of the 
enormous computational effort necessary to perform the forward wave propagation computations 
required by the algorithm. 
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However, it is likely that all of these will be overcome in time, and the straightforward 
nature of a deterministic algorithm such as this implies substantial long-term benefits to seismic 
exploration. 
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